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ABSTRACT 

Sylvia is a triple asteroid system located in the main belt. We report new adaptive optics observations of 
this system that extend the baseline of existing astrometric observations to a decade. We present the first 
fully dynamical 3-body model for this system by fitting to all available astrometric measurements. This model 
simultaneously fits for individual masses, orbits, and primary oblateness. We find that Sylvia is composed 
of a dominant central mass surrounded by two satellites orbiting at 706.5 ± 2.5 km and 1357 ± 4.0 km, 
i.e., about 5 and nearly 10 primary radii. We derive individual masses of 1 .484^[J j][4 x 10 19 kg for the primary 
(corresponding to a density of 1.29 ± 0.39 g cm" 3 ), 7.33^ 3 x 10 14 kg for the inner satellite, and9.32!|° 3 7 xlO 14 
kg for the outer satellite. The oblateness of the primary induces substantial precession and the J2 value can be 
constrained to the range of 0.0985-0. 1. The orbits of the satellites are relatively circular with eccentricities less 
than 0.04. The spin axis of the primary body and the orbital poles of both satellites are all aligned within about 
two degrees of each other, indicating a nearly coplanar configuration and suggestive of satellite formation in 
or near the equatorial plane of the primary. We also investigate the past orbital evolution of the system by 
simulating the effects of a recent passage through 3:1 mean-motion eccentricity-type resonances. In some 
scenarios this allow us to place constraints on interior structure and past eccentricities. 
Subject headings: minor planets, asteroids: general — minor planets, asteroids: individual (Sylvia) 



1. INTRODUCTION 

(87) Sylvia is a triple asteroid residing in the main belt, 
with heliocentric semi-major axis 3.5 AU, eccentricity 0.085, 
and inclination 11° relative to the ecliptic. Sylvia's outer 
satellite, named Romulus , was discovered in 2001 using the 
W. M. Keck Telescope (|Brown & Margot||200lj |Margot & 



|Brown||2001| i, and w as also detected in Hubble Space Tele- 



scope (HST) images ( Storrs et al.|200l| . The inner satellite, 
Remus, was not discovered until the advent of improved adap- 
tive optics systems in 2004 using the Eu ropean Southern Ob - 
servatory's Very Large Telescope (VLT) (Marchis et al. 2005 ). 
The diameter of the primary has been estimated at ^280 km 
throug h shape fits to adaptive optics images (Marchis et al. 
2005); this estimate is con sistent with stellar occultation ob- 
servations ( |Lin et al.|2009| l. Assuming this primary size, ap- 
proximate sizes for the individual satellites have been esti- 
mated by adopting the same albedo as the primary and mea- 
suring each satellite's brightness relative to the primary. The 
diameter estima tes are ^7 km for Remus and ^18 km for Ro- 
mulus ( |Marchis et al.|2005] >. 

Sylvia was the first triple asteroid system announced, 
even though the triple nature of 2002 CE26 was being ac- 
tively discussed during th e acquisition of the Sylvia observa- 
tions ( |Shepard et al.|[2006| l. Additional discoveries of multi- 
ples in the Solar System have followe d. They include n ear- 
Earth triples (15359 1) 2001 SN263 (prolan et al.||2008] > and 
(136617) 1 994 CC flBrozovic et al.||2009i m ain beTTtriples 
Kleopatra ( Descamps et al. ||201 1| >7 E ugenia (|Merline et al.| 
19991 IMarchis et al.||2007| r^ala m (|Merline et al.||2002[ 



Marchis et al.|2008| l,"and Minerva ( |Marchis et al.|201 1[>, and 



trans-Neptunian systems (47171) 1999 TC36 (Margot et al. 
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[2005| |Benecchi et aIJ|2lH0] i, Hau mea (|Brown et al.|[2005] 
[2006b, and the Pluto/Charon system ( [Weaver et al.|2006| >. 

Following these discoveries, characterization of multiple 
systems have unearthed a wealth of information about their 
fundamental physical properties such as masses and densi- 
ties, dynamical processes, and constraints on formation and 
evolutionary mechanisms. Such research has been possible 
because we can derive the masses of the individual compo- 
nents of a triple or higher-multiplicity system by analyzing 
their mutual gravitational interactions, which is possible in 
binary systems only when reflex motion is detecte d ( |Mar-| 
|got et al.|[2002| |Ostro et aL[ 2006; N aidu et al.||20TT] >. These 
masses in conjunction with size estimat es can provide den- 
sities. Using this method, Fang et al. ( |2011| l performed a 
detailed analysis of 2001 SN263 and 199TCC, including 
masses, densities, and dynamical evolution. Similarly, work 
on the Pluto/Charon system and dwarf planet Haumea and its 
satellites have yielded information about their physical prop- 
erties, tidal interactions, and evolutionary processes (|Lee &| 
|PeaTel[2006l |Tholen et al.||2008] |Ragozzine_ & BrownJ 2009 ). 
The high scientific return from studies of binaries and triples 
has be en reviewed by Merline et al. (2002) and Noll et al. 
( |2008l >. 

To date, no such dynamical orbit solution nor detaile d anal- 
ysis h a s been performed for Sylvia. Previous work by [Marchis| 



et al. ( 2005} approximated the actual orbits of Remus and 
Romulus with individual two-body fits that included primary 
oblateness. However, drawbacks of such methods include the 
failure to account for third-body perturbations as well as the 
inability to solve for individual component masses. Addi- 
tional resea rchers based their studies on the published two- 
body orbits (Marchis et al. 2005)> plus unspecified component 



mass assumptions to study Sylvia's long-term evolution ( Win- 
|ter et al.|20 09 ; Frouard & Comper e|2012[ ), even though com- 
ponent masses are undetermined and can span several orders 
of magnitude. 

In this work, we report additional Keck and VLT imaging 



data for Sylvia (Section |2j. Using primary-satellite separa- 
tions measured from these data plus published astrometry, we 
present a fully dynamical 3-body orbital and mass solution 
for Sylvia, by accounting for mutually interacting orbits as 
well as the primary's non-sphericity (Section B). Although 
the orbital periods of the satellites are near a 8:3 ratio, we do 
not find that the system is currently in such a resonance (Sec- 
tion |4jl. We also analyze Sylvia's short-term and long-term 
future evolution (Section BJ. Lastly, we investigate the past 
orbital evolution of Remus and Romulus by modeling pas- 
sage through the 3:1 mean-motion resonance (Section [5J. A 
summary of main conclusions is given in Section|7] 

2. OBSERVATIONS 

We report new observations of Sylvia in 2011 from both 
Keck and VLT, as well as summarize existing observations 
taken in 2001-2004 using Keck, HST, and VLT. Astrometry 
derived from these datasets are used for orbit fits described in 
the next section (Section[3]). 

2.1. New Data in 2011 

Our observations in 2011 are summarized in Table Q] In to- 
tal, we obtained 7 partial nights of service (or "queue") mode 
observing at VLT and 4 partial nights of visitor mode at Keck. 
At the VLT, we used the NA CO (NAOS-CONICA) high- 
resolution I R imaging camera (Rousset et al. 2003] |Lenzen| 
|et al.||20"03"] ) in its S13 mode using the H filter, with aplate 
scale value reported as 0.013221 arcseconds per pixefT We 
used four offset positions in a box pattern, with an integration 
time of 120 seconds per offset position. These four offset po- 
sitions sampled the four quadrants of the CCD to calibrate and 
mitigate against detector defects. At Keck, we used NIRC2 
(Near Infrared Camera 2) imaging in the H filter, with a plate 
scale of 0.009942 arcseconds per pixerl Our Keck observa- 
tions used four offset positions in a similar box pattern as with 
the VLT exposures, with an exposure time of 60 seconds per 
offset position. All observations using both VLT and Keck 
were performed with natural guide star adaptive optics, using 
Sylvia as the guide star (its apparent magnitude varied from 
V~l 1.7 to V~12.7 throughout the period of our 201 1 obser- 
vations). 

We performed basic data reduction analysis for VLT and 
Keck images. Each frame (at each dither position) is flat- 
fielded and its bad pixels are corrected using a bad pixel mask 
obtained from outlier pixels present in all frames. Sky sub- 
traction is performed by subtracting frames in pairs, where 
each frame is subtracted by another frame where the target 
has been offset. We performed subpixel 2-D Gaussian fitting 
to obtain precise centroids of the target in each sky-subtracted 
frame, and these centroids were used to align and combine all 
frames into one composite image. See Figure [1] for an exam- 
ple of a composite image obtained using Keck. 

At each observation epoch, we detected either one or both 
satellites (see Table[T]i. As the satellites orbit the primary, they 
are occasionally obscured by the bright primary and this oc- 
curs most often for the inner satellite Remus. In cases where 
only one satellite is detected, we determined the identification 
of the satellite through orbit fitting. An incorrect identifica- 
tion of a satellite is easily shown as an obvious outlier in orbit 
fits. In all cases examined here, when only one satellite is 
detected, it is the outer satellite Romulus. 

1 http://www.eso.org/sci/facilities/paranal/instruments/naco/doc 

2 http://www2.keck.hawaii.edu/inst/nirc2/genspecs.html 



For all satellite detections, we measure the astrometric po- 
sitions of the satellite relative to the primary by taking the dif- 
ference between the centroids of the primary and the satellite. 
Specifically, we measure the position angle (degrees East of 
North) and separation of the satellite relative to the primary. 
These measurements are performed on the reduced composite 
image obtained at each observation epoch. These measure- 
ments are provided in Table [2] 



Fig. 1 . — Keck H-band adaptive optics image on 201 1 Decem- 
ber 15 (corresponding modified Julian Date 55910.2129) of 
Sylvia with inner satellite Remus and outer satellite Romulus. 
In this image, the primary-Remus separation is about 0.34 
arcseconds, and the primary-Romulus separation is about 
0.57 arcseconds. 



2.2. Existing Data From 2001-2004 

Publicly-available astrometry da tasets for Sylvia include 
Keck data in 2001 ([Margot : & Brown|200T) , H ST data in 2001 
UStorrs et al.||2001) , and VLT data in 2004 ( |Marchis et"aL 
120051. For the 2004 VLT data, the paper by |Marchis et~aT 
(2005) contains their astrometric measurements of the satel- 
lites relative to the primary expressed as x—y pairs, where they 
define x and y as positive in the East and North directions, re- 
spectively. However, they failed to indicate the signs (positive 
or negative) for their x and y measurements of Remus. In 
addition, at one epoch (MJD 53253.1738) their published as- 
trometry give Remus an implausibly large separation from the 
primary. To fix these inaccuracies, we fit orbits to the astro- 
metric points and have determined the correct signs for their 
measurements (note that we define x to be positive in the West 
direction). At epoch MJD 53253.1738, it appears that they 
have confused Remus and Romulus, and so we have swapped 
measurements for these two bodies at that epoch. These cor- 
rections, along with astrometry for the Keck and HST data in 
2001, are given in Table [3] 

In total, our baseline of observations spans a decade from 
2001 to 2011. For Remus, which is never visible in data ob- 
tained prior to 2004, our 201 1 observations add an additional 
8 epochs to the existing 12 epochs for a total of 20 epochs of 
astrometry, extending the baseline of observations from about 
1 month to 7 years. For Romulus, our 201 1 observations add 
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TABLE 1 
Summary of 201 1 Observations 



UT Date 


MJD 


Filter 


Telescope 


Detections 


2011 Oct 7 


55841.1097 


H 


VLT 


Remus, Romulus 


2011 Nov 6 


55871.0856 


H 


VLT 


Romulus 


2011 Nov 8 


55873.1289 


H 


VLT 


Remus, Romulus 


2011 Nov 10 


55875.0451 


H 


VLT 


Remus, Romulus 


2011 Nov 15 


55880.0256 


H 


VLT 


Romulus 


2011 Nov 16 


55881.0466 


H 


VLT 


Remus, Romulus 


2011 Nov 20 


55885.0460 


H 


VLT 


Romulus 


2011 Dec 15 


55910.2129 


H 


Keck 


Remus, Romulus 


2011 Dec 15 


55910.2288 


H 


Keck 


Remus, Romulus 


2011 Dec 15 


55910.2698 


H 


Keck 


Remus, Romulus 


2011 Dec 16 


55911.1877 


H 


Keck 


Romulus 


2011 Dec 16 


55911.2510 


H 


Keck 


Romulus 


2011 Dec 17 


55912.2615 


H 


Keck 


Remus, Romulus 


2011 Dec 18 


55913.1972 


H 


Keck 


Romulus 


2011 Dec 18 


55913.2035 


H 


Keck 


Romulus 



Summary of our 2011 adaptive optics observations at VLT and Keck. 
Epochs are provided in Universal Time (UT) dates as well as the Modi- 
fied Julian Date (MJD). Remus is the inner satellite and Romulus is the 
outer satellite. 



TABLE 2 
ASTROMETRY FROM 201 1 DATA 



Satellite 


MJD 


PA 


Sep. 


X 


y 








(deg) 


(arcsec) 


(arcsec) 


(arcsec) 


(arcsec) 


Remus 


55841.1097 


265.07 


0.392 


0.391 


-0.034 


0.0132 


Remus 


55873.1289 


75.25 


0.226 


-0.218 


0.057 


0.0132 


Remus 


55875.0451 


270.33 


0.342 


0.342 


0.002 


0.0132 


Remus 


55881.0466 


87.23 


0.381 


-0.381 


0.018 


0.0132 


Remus 


55910.2129 


268.16 


0.341 


0.341 


-0.011 


0.0099 


Remus 


55910.2288 


267.49 


0.343 


0.342 


-0.015 


0.0099 


Remus 


55910.2698 


267.80 


0.325 


0.325 


-0.012 


0.0099 


Remus 


55912.2615 


84.22 


0.349 


-0.348 


0.035 


0.0099 


Romulus 


55841.1097 


264.09 


0.702 


0.698 


-0.072 


0.0132 


Romulus 


55871.0856 


92.92 


0.369 


-0.368 


-0.019 


0.0132 


Romulus 


55873.1289 


271.41 


0.606 


0.606 


0.015 


0.0132 


Romulus 


55875.0451 


88.92 


0.643 


-0.643 


0.012 


0.0132 


Romulus 


55880.0256 


284.49 


0.183 


0.177 


0.046 


0.0132 


Romulus 


55881.0466 


266.17 


0.671 


0.670 


-0.045 


0.0132 


Romulus 


55885.0460 


264.34 


0.357 


0.355 


-0.035 


0.0132 


Romulus 


55910.2129 


264.30 


0.571 


0.568 


-0.057 


0.0099 


Romulus 


55910.2288 


263.66 


0.572 


0.568 


-0.063 


0.0099 


Romulus 


55910.2698 


265.11 


0.541 


0.539 


-0.046 


0.0099 


Romulus 


55911.1877 


95.41 


0.392 


-0.391 


-0.037 


0.0099 


Romulus 


55911.2510 


92.09 


0.446 


-0.446 


-0.016 


0.0099 


Romulus 


55912.2615 


78.76 


0.423 


-0.415 


0.082 


0.0099 


Romulus 


55913.1972 


272.80 


0.528 


0.528 


0.026 


0.0099 


Romulus 


55913.2035 


272.09 


0.521 


0.521 


0.019 


0.0099 



Astrometry measured from 2011 data for Remus (inner) and Romulus (outer) 
at specific MJD epochs. We measured the position angle (PA; degrees East of 
North) and the separation of each satellite relative to the primary. These values 
are converted to positions x and y, where positive x direction is towards the West, 
and positive y direction is towards the North. We assign the instrument plate scale 
as the positional uncertainty a for x and y. 
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an additional 15 epochs to the existing 30 epochs for a total of 
45 epochs of astrometric measurements, extending the base- 
line of observations from 3 years to 10 years. All of these 
astrometry positions, given in Tables |2l and [3] are used for 
dynamical three-body orbit fits described in the next section. 

3. ORBITAL AND MASS SOLUTION 

We fit a fully dynamical 3 -body model to the astrometric 
measurements described in the previous section, taking into 
account mutually interacting orbits. Our model simultane- 
ously fits for 16 parameters, including a set of 6 orbital param- 
eters per satellite (semi-major axis, eccentricity, inclination, 
argument of pericenter, longitude of the ascending node, and 
mean anomaly at epoch), 3 masses for the three bodies, and 
the parameter representing the oblateness of the primary. 
The semi-major axis and eccentricity describe the size and 
shape of the orbit, both the inclination and longitude of the 
ascending node describe the orientation of the orbit, the argu- 
ment of pericenter describes the location of pericenter (min- 
imum radial distance of the orbit), and the mean anomaly at 
epoch can be used to determine the location of the satellite in 
its orbit at a particular time. 

We describe these fitted parameters in more detail. The 
orbital elements of each satellite are relative to the primary 
body, and are defined with respect to the Earth equatorial ref- 
erence frame of epoch J2000. Given that these orbital ele- 
ments change over time due to perturbations in the three-body 
system, these are defined as osculating orbital elements. They 
are valid at a specific epoch, MJD 53227.0, corresponding to 
UT date 2004 August 10 00:00. The fitted masses are de- 
rived assuming G = 6.67 x 10~ n m 3 kg -1 s~ 2 for the gravi- 
tational constant. The primary's non-spherical nature, which 
can introduce additional non-Keplerian effects, is represented 
by an oblateness coefficient J^- The distribution of mass 
within the primary can be represented by terms in a spheri- 
cal harmonic expansion of its gravitational potential, and the 
quadmpole term J2 is the lowest-order gravitational moment. 
J2 is related to the primary's three principal moments of iner- 
tia (C > B > A) as 



h = 



C- l -{A+B) 
MR 2 ■ 



(1) 



where the denominator is a normalization fact or including the 
primary's mass M and equatorial radius R (i.e., Murray & Der- 
|mott|199 9). In all of our fits, the radius R is assumed to be 140 
km. The inclusion of J2 in our fits also requires a primary spin 
pole direction to be specified, and typically we fix the primary 
pole to the orbit pole of the most massive satellite. 

With a total of 16 parameters and 130 data measurements 
(Tables [2] and [3J, we have 1 14 degrees of freedom (number 
of data points minus the number of parameters). We adopt 
a least-squares approach to this problem by minimizing the 

chi-square x 2 statistic, where \ 2 = J^(<3; _ C,) 2 /of for a set 

i 

of N(l < i < N) observations with <r, uncertainties and ob- 
served O, and computed C, values. We utilize a Levenberg- 
Marquardt non- linear least- squares algorithm written in IDL 
called mpf it (Markwardt 2009). With 16 parameters, this 
is a very computationally intensive problem, given the large 
amount of parameter space to explore as well as the compu- 
tationally expensive 3 -body orbital integrations that need to 
be performed. Especially for least-squares problems with a 



large number of parameters, it is impossible to guarantee that 
a global \ 2 minimum has been found. More often than not, 
the minimization procedure converged on a local minimum 
and therefore it was necessary to re-fit with different starting 
conditions. In total, we started the fitting procedure with tens 
of thousands of sets of starting conditions, and we performed 
up to 20 iterations (equaling hundreds of x 2 evaluations) for 
each set of starting conditions. 

These starting conditions included plausible ranges of pa- 
rameter space for fitted parameters. We explored all possible 
values of eccentricity (0-1), orbital angles (0-360°), semi- 
major axes (500-900 km for Remus and 1100-1500 km for 
Romulus), and J2 values (0-0.2). Starting values ranged from 
on the order of 10 18 to 10 20 kg for the primary's mass and 
from on the order of 10 13 to 10 17 kg for each satellite's mass. 
Ranges for satellite masses covered all possible mass values 
by sampling various size (Remus: ~5-9 km in di ameter, Ro- 
mulus: ^14-22 km in diameter; Marchis et al. 2005) and 
density (0.1-10 g cm" 3 ) ranges. 

In addition, we also explored various primary spin axis ori- 
entations for the primary. This was possible with this data set 
because the perturbations due to the oblateness of the primary 
are detectable, and because those perturbations depend on the 
spin axis orientation. These effects are captured by three fitted 
parameters: two for the primary spin axis orientation, and one 
for the value of J 2- We systematically explored the entire ce- 
lestial sphere for the primary spin axis orientation, but we also 
tested specific poles that had been favored by previous stud- 
ies. These specific spin axis orientations include RA=355° 
and DEC=82° (close to s atellite orbital poles) suggested from 
light curve analysis (Kaasalai nen et al.||2002| l , RA=68° and 
DEC=78° f rom a compilation of previous data ( |Kryszczynska| 
et al.|2007) , and RA=100 ° and DEC=62° derived using adap- 



tive optics imaging data (Dm mmond & Christou|2008). From 
these fits, we find that primary spin poles misaligned with 
satellite orbit poles do not provide go od solutions (such as 
poles bylRr yszczyns ka et aL"1 ( 2007| ) and Drummond & Chris- 



tou| ( |2008| l). Instead, we determined that the best-fit spin axis 
direction was nearly aligned with the satellites' orbital poles, 
which are almost coplanar. As a result, for nearly all fits, we 
aligned the primary's spin pole to the orbital pole of the most 
massive satellite. The fact that both satellites orbit in or near 
the equator of the primary provides an important constraint on 
satellite formation mechanisms. 

For each set of starting conditions, our orbit-fitting method 
proceeded as follows. First, we performed N-body numerical 



integrations using the Mercury integration package ( [Cham- 
bers) 19 99), which takes into account mutually interacting or- 



bits as well as the effects due to primary oblateness. We 
used a Bulirsch-Stoer algorithm for our integration method, 
which is computationally slow but accurate, and we chose 
an initial time step that samples finer than l/25th of the in- 
nermost orbital period. These 3-body integrations need to 
cover all epochs of observation, which span about a decade 
(2001-2011). From these integrations, we determined the 
positions and velocities of each satellite relative to the pri- 
mary at all epochs of observation (corrected for light travel 
time) by interpolation. The length and resolution of these 
integrations were the limiting factors in the computational 
speed of each minimization. Second, we obtained the vec- 
tor orientation of Sylvia's position relative to an observer on 
Earth for all epochs of observation (again, corrected for light 
time), taking into account aspect variations due to geocen- 
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TABLE 3 

Existing Astrometry From 2001-2004 Data 



Satellite 


MJD 


PA 


Sep. 


X 


y 


(7 


Reference 






(deg) 


(arcsec) 


(arcsec) 


(arcsec) 


(arcsec) 




Remus 


53227.3042 






-0.411 


0.002 


0.0132 


3 


Remus 


53249.2470 






-0.239 


-0.101 


0.0132 


3 


Remus 


53249.2532 






-0.228 


-0.101 


0.0132 


3 


Remus 


53251.2985 






0.246 


0.107 


0.0132 


3 


Remus 


53252.3627 






0.421 


0.000 


0.0132 


3 


Remus 


53253.1738 






-0.445 


-0.025 


0.0132 


3 


Remus 


53255.1091 






0.435 


0.002 


0.0132 


3 


Remus 


53256.2886 






0.268 


-0.072 


0.0132 


3 


Remus 


53261.1432 






-0.394 


0.033 


0.0132 


3 


Remus 


53261.2298 






-0.430 


-0.008 


0.0132 


3 


Remus 


53263.2146 






0.412 


0.004 


0.0132 


3 


Remus 


53263.2202 






0.416 


0.001 


0.0132 


3 


Romulus 


51958.4810 


97.00 


0.564 


-0.559 


-0.069 


0.0168 


1 


Romulus 


51959.3660 


60.30 


0.425 


-0.369 


0.211 


0.0168 


1 


Romulus 


51959.4200 


54.10 


0.383 


-0.310 


0.225 


0.0168 


1 


Romulus 


51960.4010 


271.90 


0.615 


0.614 


0.020 


0.0168 


1 


Romulus 


51962.4070 


88.30 


0.696 


-0.696 


0.021 


0.0168 


1 


Romulus 


51963.5700 


306.00 


0.330 


0.267 


0.194 


0.0250 


2 


Romulus 


53227.3042 






-0.377 


0.144 


0.0132 


3 


Romulus 


53246.3105 






-0.785 


-0.097 


0.0132 


3 


Romulus 


53246.3659 






-0.755 


-0.117 


0.0132 


3 


Romulus 


53249.2470 






-0.547 


0.139 


0.0132 


3 


Romulus 


53249.2532 






-0.555 


0.136 


0.0132 


3 


Romulus 


53249.3516 






-0.654 


0.109 


0.0132 


3 


Romulus 


53251.2985 






0.763 


-0.058 


0.0132 


3 


Romulus 


53252.3627 






0.156 


0.214 


0.0132 


3 


Romulus 


53253.1738 






-0.791 


0.049 


0.0132 


3 


Romulus 


53253.3445 






-0.834 


-0.018 


0.0132 


3 


Romulus 


53254.1603 






-0.172 


-0.214 


0.0132 


3 


Romulus 


53255.1091 






0.835 


0.003 


0.0132 


3 


Romulus 


53255.3928 






0.793 


0.105 


0.0132 


3 


Romulus 








O 979 


O 909 


o ni 39 

U.UljZ 


~i 


Romulus 


53259.2030 






0.683 


0.165 


0.0132 


3 


Romulus 


53261.1432 






-0.546 


-0.186 


0.0132 


3 


Romulus 


53261.2298 






-0.449 


-0.205 


0.0132 


3 


Romulus 


53262.1602 






0.724 


-0.073 


0.0132 


3 


Romulus 


53262.2759 






0.786 


-0.035 


0.0132 


3 


Romulus 


53262.2815 






0.791 


-0.032 


0.0132 


3 


Romulus 


53263.2146 






0.221 


0.230 


0.0132 


3 


Romulus 


53263.2202 






0.215 


0.231 


0.0132 


3 


Romulus 


53297.0193 






-0.752 


-0.038 


0.0132 


3 


Romulus 


53298.0064 






0.101 


-0.186 


0.0132 


3 



Existing astrometry measured from 2001-2004 data for Remus (inner) and Romulus (outer) 
at specific MJD epochs taken from [1] Margot & Brown 12001 1, [2] Storrs et al. |200U, and 
[3] Marchis et al. (2005) (the latter with corrections; see text). An ellipsis (...) means that the 
value was not reported. Measurements include the position angle (PA; degrees Earth of North) 
and the separation of each satellite relative to the primary. These values can be converted 
to positions x and y, where positive x direction is towards the West, and positive y direction 
is towards the North. As in Table [2] we assign the instrument plate scale as the positional 
uncertainty a for x and y. 



5 



trie distance variations and Sylvia's motion across the sky. 
Third, we used these orientations to project and compute 
primary-satellite separations on the plane of the sky at each 
observation epoch. These computed separations were com- 
pared with our observed separations (Tables [2]and[3]) to deter- 
mine the x 2 goodness-of-fit statistic. These methods benefit 
from the heritage of over a decade of work on orbit fitting 
as well as our work on fitting 3 -body models to datas ets of 
triples, including near-Earth asteroids (Fang et al.|201 1] >. 

Table |4] shows the best- fit orbit solution. The cni-square is 
73.55, and with 1 14 degrees of freedom, this corresponds to 
a reduced chi-square of 0.6452. This indicates that the fit is 
very good, and that uncertainties were likely slightly overes- 
timated. This best-fit solution is also visually illustrated in 
an orbit diagram in Figure [2] where the orbits are projected 
onto the primary's equatorial plane. Residuals of the best-fit 
solution are shown in Figure [3] for Remus and Figure H] for 
Romulus, where each figure's panel shows the residuals for 
a particular year of observation (2001, 2004, or 2011). The 
residual is defined as Xi = (O,-— Q)/<7j for observation epoch ;. 

There are two types of la uncertainties given in Table Hj 
formal and adopted. The formal la uncertainties are orv 
tained using the covariance matrix resulting from the least- 
squares fitting procedure. These formal la uncertainties may 
not always be accurate representations of the actual uncertain- 
ties. Accordingly, the adopted la uncertainties are obtained 
through a more rigorous method by determin ing each parame - 
ter's la confidence levels (e.g., |Cash|1976| [Press et al.|1992| l. 
To determine each parameter's uncertainties, we hold the pa- 
rameter fixed at a range of plausible values while simultane- 
ously fitting for all other parameters. Since one parameter 
is held fixed at a time, a la confidence region is prescribed 
by the range of solutions that yield chi-square values within 
1.0 of the lowest chi-square. This is a computationally in- 
tensive process, and we performed this method to determine 
uncertainties for the primary's pole (R.A. and Dec.) as well as 
for each fitted parameter with the exception of the arguments 
of pericenter and mean anomalies at epoch. We consider the 
adopted uncertainties to be more accurate representations of 
the actual uncertainties than the formal uncertainties. We note 
that these adopted uncertainties do not preclude any system- 
atic errors that may have occurred, such as during the mea- 
surement of astro me try from the images. 

Most solve-for parameters are well constrained by the data, 
with formal uncertainties of 10% or less, and as low as ~1% 
for the mass and oblateness of the primary. One notable ex- 
ception is the mass of Romulus. The formal uncertainty on 
this parameter amounts to ^60% of the nominal value, in- 
dicating that it is not well constrained by the data. As for 
orbit pole orientations, they are determined to ~1.5 degrees 
uncertainty (la) for Remus and ~1 degree for Romulus, and 
the primary spin axis orientation is determined to ~1 degree. 
Given the near-circular nature of the orbits, the arguments of 
pericenter lu and hence mean anomalies at epoch M are not 
well-defined (but the satellite positions, or w+M, are in fact 
well-defined). 

The best-fit orbital solution indicates that the satellites fol- 
low relatively circular orbits at semi-major axes of about 5 
and nearly 10 primary radii. The mutual inclination (rela- 
tive inclination) between the orbital planes of the satellites is 
0.56°^o5. The low mutual inclination indicates that the orbital 
planes are nearly coplanar. This alignment supports a satellite 
formation mechanism in which the satellites would tend to re- 



TABLE 4 

Best-fit Parameters and la Uncertainties 







Best-fit 


Formal la 


Adopted la 






Remus (inner): 




IVldSS \1V 


kg) 


7.333 


±0.7172 


+4.667 
-2.333 


ci (km) 




706.5 


±0.007231 


+2.488 
-2.512 


c 




0.02721 


±0.009962 


+0.01279 
-0.01221 


i (deg) 




7.824 


±0.6665 


+0.6763 
-0.8237 


(jj ( (\pa\ 




357.0 


±15.14 




\L [peg) 




94.80 


±5.000 


+5.198 
-5.802 


M f rWi 
ivi yucgj 




261.0 


±13.43 




r (clays } 




1.373 




+0.01036 
-0.009771 






Romulus (outer): 




Mass HO 14 


kg) 


9.319 


±5.406 


+20.68 
-8.319 


a (km) 




1357 


±0.05918 


+3.984 
-4.016 






0.005566 


±0.004268 


+0.005434 
-0.003566 


, 

i (aegj 




8.293 


±0.2099 


+0.2069 
-0.2931 


uj (deg) 




61.06 


±18.74 




H (deg) 




92.60 


±1.339 


+2.903 
-1.597 


M(deg) 




197.0 


±18.75 




P (days) 




3.654 




+0.02544 
-0.02366 






Primary: 




Mass (10 19 


kg) 


1.484 


±0.0001659 


+0.01601 
-0.01399 


h 




0.09959 


±0.0008384 


+0.0004148 
-0.001085 


R.A. (deg) 




2.597 


±1.339 


+3.403 
-1.597 


Dec. (deg) 




81.71 


±0.2099 


+0.2931 
-0.7069 



Best-fit parameters including individual masses, orbital pa- 
rameters (semi-major axis a, eccentricity e, inclination i, ar- 
gument of pericenter u, longitude of the ascending node f!, 
mean anomaly at epoch M), primary oblateness J2, and pri- 
mary spin pole (R.A. and Dec). These orbital elements are 
valid at epoch MJD 53227 in the equatorial frame of J2000. 
We derived an effective orbital period P from the best-fit val- 
ues of semi-major axis and mass of the considered satellite 
plus all interior masses. Two types of uncertainties are listed: 
formal la statistical errors are derived from the least-squares 
covariance matrix, and adopted la errors are obtained for se- 
lect parameters through a more rigorous method (see text). 
For parameters lu and M with no adopted la errors, we rec- 
ommend using the formal la errors. 

main close to the equatorial plane of the primary (e.g., a sub- 
catastrophic collision). The alignment is likely indicative of 
formation conditions rather than evolution, as tidal damping 
of inclinations can be a lengthy process. Assuming various 
models (monolithic or rubble pile) for tidal dissipation, inner 
satellite Remus could take on the order of 10 8 years up to the 
age of the Solar System to damp from 2 degrees to 1 degree. 

Information about size, shape, and density can be derived 
from our best-fit orbital solution. We obtain a density of 
1.29±0.39 g cm -3 for the primary by assuming a diameter 
of 280 km (lacking realistic error bars on the size of the pri- 
mary, we assumed volume uncertainties of 30%). The den- 
sity error is dominated by the volume error since the mass of 
the primary is known to ~1%. The primary is also oblate, 
with a well-constrained J2 value of about 0.09959 which cor- 
responds to an axial ratio c/a = 0.7086 if we assume equato- 
rial symmetry and uniform density. Size estimates for Remus 
and Romulus can be obtained by assuming that they have a 
bulk density equal to that of the primary, and by considering 
the adopted la confidence interval of satellite masses. We 
find radii of ~4.5— 6.1 km for Remus and ^2.6-8.2 km for 
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Fig. 2. — Diagram of best-fit orbits for satellites Remus (in- 
ner) and Romulus (outer), projected onto the primary's equa- 
torial plane. These orbits show the actual trajectories from 
numerical integrations. The relative sizes of the bodies are 
shown to scale using green circles, assuming these diameters: 
10.6 km for Remus, 10.8 km for Romulus, and 280 km for 
the primary. All bodies are located at their positions at MJD 
53227 with the primary centered on the origin. 
Romulus. These ranges would have to be modified if the den- 
sity of the primary or of the satellites was different from the 
nominal value assumed here, but only slightly as the depen- 
dence is p" 1 / 3 . 

We compare our bes t- fit solution in Table |4| to the solution 
previously reported by Marchis et al.| (|2"005). We find close 
agreement in semi-major axes (within uncertainties). The ec- 
centricities are marginally consistent. Orbital plane orienta- 
tions differ by about <2 degrees. The largest discrepancy be- 
tween our orbital solutions is the value of ]%. Our fits yield a 
very well-const rained J2 va lue with a ler confidence range of 
0.0985-0.1. Marchi s~etaLl ( |2005| > report two estimates for J 2 
of 0. 17 ± 0.05 and 0.18 ± 0.01, inconsistent with our value- 
Using axial ratios from lightcurve analysis (Kaasalaine n et al.| 
2002 ) and a uniform density assumption, we find J 2 = 0.1, in 
excellent agreement with our dynamical value. 

There are several possibilities to explain th e discrepancies 
betwee n our orbital solutions and those of Marchi s et al.l 
( |2005[ ). First, our orbital fits are based on a much longer 
baseline of observations (2001-201 1) than their dataset (only 
2004). Second, our orbital solution is the result of fully dy- 
namical, N-body orbital fits that simultaneously fitted for all 
parameters in the system, using numerical integrations taking 
into account mutually int eracting orbits and primary oblate- 
ness. The fit obtained by |Marchis et al.| ( [2005| l used two-body 
approximations (one satellite's orbit is fit at a time, ignoring 
effects by the other satellite). These differences in dataset 
and technique allowed us to obtain better-constrained orbital 
parameters with smaller uncertainties as well as individual 
masses, which were not previously known. 

4. EXAMINATION OF MEAN-MOTION RESONANCE OCCUPANCY 

The orbital periods (ratio ss 2.661) of our best-fit solution 
in Table [4] have a ratio near 8:3 (ratio « 2.667). To deter- 
mine resonance occupation, we search for librating resonance 
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Fig. 3. — Residuals for separations in the West and North di- 
rections corresponding to the best-fit orbit (Table Q for inner 
satellite Remus. The two panels represent distinct epochs of 
observations (2004, 2011). 

arguments using a genera l form of the resonance argument 
( |Murray & Dermott|1999| l 

= j\ A2 + 72A1 + ; 3 a3 2 + + ;' 5 ^2 + 76^1 ■ (2) 

In Equation p]), <\> is the resonant argument or angle, A is the 
mean longitude, 03 is the longitude of pericenter, and Q is the 
longitude of the ascending node. Subscripts 1 and 2 repre- 
sent the inner and outer satellites, respectively. The ji val- 
ues (where i = 1-6) are integers and their sum must equal 
zero (d' Alembert's rule). For the fifth-order 8:3 mean-motion 
resonance, j\ = -8 and ]2 = 3 so we search through integer 
values (-30 to +30) of the remaining _/',- values to determine 
if there is libration of the resonant argument over timescales 
ranging from 10 to 100 years. To perform this search, we de- 
termined the evolution of the relevant angles in Equation d2]i 
using 3-body numerical integrations with the best-fit orbital, 
mass, and J2 solution in Table [4] We do not find any librat- 
ing resonance arguments, and therefore we conclude that the 
current system is not in the 8:3 mean-motion resonance. 

5. SHORT-TERM AND LONG-TERM STABILITY 

In this section, we discuss the results of forward N-body in- 
tegrations of the best-fit orbital solution at MJD 53227 given 
in Table |4] We perform short-term (50 yr) and long-term 
(1 Myr) simulations to determine how the orbital elements 
fluctuate with time and to assess the stability of this three- 
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Fig. 4. — Residuals for separations in the West and North di- 
rections corresponding to the best-fit orbit (Table [4]l for outer 
satellite Romulus. The three panels represent distinct epochs 
of observations (2001, 2004, 201 1). 

body system. Both short-term and long-term integrations 
are performed u sing a Bulirsch-Stoer algorithm in Mercury 
( ] Cham bers 1999) with an initial timestep of 0.05 days, and 
include the gravitational effects of the three bodies and the 
primary's oblateness. For long-term integrations, we also in- 
clude the effect of solar gravitational perturbations. 

The results of our short-term integrations are shown in Fig- 
ure [5] This figure illustrates how the semi-major axes and 
eccentricities of both satellites evolve over a span of 50 years. 



From this figure, we can compute the mean value of orbital el- 
ements, in contrast to the osculating orbital elements provided 
in Table[4]that are valid at the specific epoch of MJD 53227.0. 
The semi-major axes for both Remus and Romulus have small 
oscillations spanning less than 1 km, and have mean values of 
706.57 km and 1356.83 km, respectively. The mean eccen- 
tricity values are 0.029 for Remus and 0.0074 for Romulus. 
The eccentricity variations are especially apparent for Remus, 
whose eccentricity can vary from 0.023 to 0.035. These short- 
period fluctuations are due to the effect of the oblateness J2 of 
the primary. The force due to the primary's gravitational field 
can be modified to account for primary J 2, and this modified 
force affects a satellite's orbit by inducing short-period fluctu- 
ations in the semi-major axis, mean motion, eccentricity, and 
mean anomaly. Its effec t on eccentricity can be mathemati- 
cally approximated as ( |Brouwer|1959| |Greenberg|1981[ ) 



Ae a 3J 2 (R p /a) 2 



(3) 



where Ae is the maximum eccentricity excursion from min- 
ima to maxima and R p is the primary's radius. Plugging in 
values of J 2 = 0.09959, R p = 140 km, and a = 706.57 km (Re- 
mus) and 1356.83 km (Romulus), we compute Ae w 0.0117 
for Remus and Ae m 0.00318 for Romulus. These values are 
consistent with the maximum excursions seen in numerical 
simulations that are plotted in Figure [5] Since Remus has a 
smaller separation (~5 R p ) from the primary than Romulus 
(~9.7 R p ), its perturbation by primary J 2 is stronger, hence 
the larger eccentricity variations seen in Figure [5] As for pre- 
cession of the orbital planes, there is significant precession 
due to J 2 and the presence of the other satellite. For example, 
Remus' longitude of pericenter precesses ^560° per year due 
to J 2 and ~1° per year due to Romulus. 

The results from our long-term integrations show that this 
triple system is stable over 1 Myr, and the variations in semi- 
major axis and eccentricity do not noticeably exceed the fluc- 
tuations shown in FigureBjfor the short-term integrations. Ac- 
cordingly, we find that Sylvia is in a very stable configuration, 
as suggested by its near circular and coplanar orbital state. We 
find that the inclusion of the Sun's gravitational effect does 
not appreciably affect the semi-major axes and eccentricities 
of Sylvia's satellite orbits. If we consider Sylvia's Hill sphere 
of gravitational influence, defined as rem = Uq(M /(SM^)) 1 ' 3 
(a Q is the heliocentric semi-major axis, M is the primary's 
mass, and M Q is the mass of the Sun), we calculate that the 
satellites orbit at ~1% and ^2% of the Hill radius. As a result, 
they are both well within the primary's sphere of gravitational 
dominance over the Sun. 

Our stability results are in agree ment with tw o previous in- 
vestigations on Sylvia's stability. Winter et al. (2009i per- 
formed stability analyses of the system, including the effects 
of the Sun and Jupiter. They find that Sylvia is not stable 
unless the primary has at least a minimal amount of oblate- 
ness (0.1% of their assumed primary J2 of 0.17). They show 
that the inclusion of primary oblateness gives rise to a secu- 
lar eigenfrequency that is much faster than those induced by 
other gravitational perturbations, which provides a sta bilizing 
effec t on th e satellites' orbital evolution. Frou ard & Com-I 
pere (2012) investigated Sylvia's short-term (20 years) and 
long-term evolution (6600 years) including the primary's non- 
sphericity (assuming J2 ~ 0. 14) and solar perturbations. They 
also varied the semi-major axis and eccentricity of the satel- 
lites' orbits to determine the extent of their stability zones. 
They find that the current configuration of the system lies in 
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Fig. 5. — Variation of semi-major axis and eccentricity over a 
50-year time span. 

a very stable zone. Authors from both papers ( |Winter et al.| 
|2009 |Frouard & Com pere 2012) mention that the effect of 
Jupiter is negligible compared to the effect of the Sun. 

6. EVOLUTION OF ORBITAL CONFIGURATION 

In this section, we investigate the past orbital evolution of 
Remus and Romulus. We find that tidal perturbations can 
cause the orbits to evolve and to cross mean-motion reso- 
nances. This resonance passage may perturb orbits by in- 
creasing eccentricities. First we discuss how tidal processes 
likely caused the satellites to encounter the 3:1 mean-motion 
resonance in their past, then we describe our numerical mod- 
eling methods, and lastly we present plausible past evolution- 
ary pathways as suggested by our simulation results. 

6.1. Tidal Theory 

Tidal evolution can cause the semi-major axis of an orbit 
to expand due to tides raised on the primary by its satellite 
( |Goldreich|[T963l |Goldreich & Soter|[l966l >. Tides raised on 
the satellite by the primary have an insignificant effect on 
the semi-major axis. The rate of semi-major axis evolution 
is given as 



da k p M s ( R p 
dt ~ Q p Mp 



(4) 



where a is the semi-major axis, k is the tidal Love number, Q 
is the tidal dissipation factor, M is the mass, R is the radius, 
and n is the mean motion. The subscripts p and s denote quan- 
tities for the primary and satellite, respectively. It is likely 
that tidal evolution is causing the orbits of Remus and Romu- 
lus to expand, and that their orbits were in a more compact 
configuration in the past. We discuss the relative importance 
of tides compared to another important evolutionary process 
(BYORP) at the end of this subsection. 

We expect that orbital expansion by tides is causing the rela- 
tive orbits of Remus and Romulus to slowly converge towards 
each other. Two orbits are converging if a'i/a'2 is greater than 
one, and here subscripts 1 and 2 represent Remus (inner) and 
Romulus (outer), respectively. Using Equation Q, we can 
express this criterion as 



a i 



Mi 

Mi 



11/2 
1 > 1. 



(5) 



Assuming best-fit values for the semi-major axes of Remus 
and Romulus (Table |4j, their orbits are currently converging 
as long as their mass ratio satisfies Mi /M2 > 0.0276. Taking 
into account the range of the ler adopted confidence interval 
in masses (from Table [4j, we find that this ratio is satisfied in 
all cases and therefore we expect that their orbits are converg- 
ing. The steep dependence of tidal evolution on semi-major 
axis causes the orbit of Remus to expand much faster than the 
orbit of Romulus. Given that their orbits are slowly converg- 
ing over time, we can determine the most recent mean-motion 
resonance passage encountered by the satellites. By consider- 
ing all first, second, third, and fourth order resonances (where 
a p + q:p resonance is gth order), we expect that the most re- 
cent resonance encountered by the system is the second-order 
3:1 resonance. Accordingly, in our analysis here we focus on 
the 3:1 resonance and its effect on the satellites' eccentricity 
evolution. 

We describe how tidal evolution causes the eccentricity of 
an orbit to increase or decrease. This is a competing process 
between the opposing effects of tides raised on the primary 
(eccentricity increases) and tides raised on the satellite (eccen- 
tricity decreases). These two oppo sing effects are contained 
in two term s in the equation ( Goldreich| |1963| [Goldreich & 
|Soter|1966) 



de 
dt 



57 k p M s 



8 QpM p \ a 



Rr 



21 k s M p 



(6) 



which gives the evolution of eccentricity e. The variables and 
subscripts in Equation ([6]) are the same as for Equation 

We discuss models for the tidal Love number k in Equations 
Q and |6]) used for calculations of tidal evolution in asteroids. 
These models are dependent on the asteroid's radius R. First, 
we consider the monolith model where 



1.5 



l + 2x 10 s 



1 km 
R 



(7) 



and this model is appropriate f or asteroids that are ideali zed as 
uniform bodies with no voids (Gol dreich & Sari|2009|l. Sec- 
ond, we consider a rubble pile model by Goldreich & Sari 



(2009) with the following Love number formalism: 
£«lxl0" 5 ' 



1 km 



(8) 



9 



Rubble pile models are appropriate for asteroids idealized as 
gravitational aggregates. Another rubble pile model is given 
by Jacobson & Scheeres ( 2.01 1| > where 



k^2.5x 10" 



1 km 
R 



(9) 



Equation Q was obtained by fitting to the configurations of 
known asteroid binaries and assuming they are in an equi- 
librium state with tidal and BYORP effects canceling each 
other. Comparison between these three Love number models 
are given in Figure [6] 
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Fig. 6. — Tidal Love number models (Equations UWfy as a 
function of radius. Intersections between the k oc l/R model 
and the other two models occur at 1.58 km and 14.94 km. 

Here we assume that tidal evolution is a dominant process 
and we do not analyze other evolutionary processes such as 
BYORP perturbations. BYORP is a radiative effect that is 
predicted to cause o rbital evolution of a synchronou s satellite 
on short timescales ( Cuk & Burns|200"5 Cuk|2007 1, and has 
not been observationally verified yet. BYORP effects domi- 
nate at larger semi- major axes, and tides domi nate at smaller 
semi-major axes (J acobson & Scheeres||20lT] ). Tides result 
in an increase of the semi-major axes, but BYORP can re- 
sult in an increase or a decrease, depending on the shapes 
of the satellites, which are unknown. For the purpose of 
our evolution calculations, the relevant semi-major axis rate 
is the relative migration rate of the satellites, and we com- 
pare the contributions by tides and BYORP by calculating 
|(atid es a-atides,2)/(4yoip,i-flbyoi P ,2)|, with subscripts l=inner 
and 2=outer. We find that this quantity is ~ 1.4— 3.3 (taking 
into account whether BYORP acts in the same or opposite 
directions for both satellites), and therefore we expect that 
tidal evolution will dominate the relative rate of the satel- 
lites' orbits as they converge. This calculation assumes cur- 
rent semi-major axes, Q p = 100, and a monolith tidal Love 
number model for the primary. If BYORP causes the orbits 
of the satellites to converge, then joint evolution by tides and 
BYORP will result in a higher relative migration rate than 
the tides-only evolution considered here. If BYORP works 
against tides, then their joint evolution will be slower. Given 
the order-of-magnitude uncertainties already inherent in un- 
known tidal parameters such as Q and k as well as uncertain- 
ties in whether BYORP will expand or contract the satellites' 
orbits, we do not include the effects of BYORP in our simu- 
lations. 



6.2. 3:1 Eccentricity-type Resonances 

The 3:1 mean-motion resonance is the most recent low- 
order resonance encountered by Remus and Romulus, and 
here we briefly describe the relevant eccentricity-type res- 
onances that can affect orbital eccentricities. We do not 
consider 3:1 inclination-type resonances in our simulations. 
There are 3 eccentricity-type resonances for the 3:1 mean- 
motion resonance: e\ resonance which perturbs only the outer 
satellite's eccentricity, <?ie 2 mixed resonance which perturbs 
both satellites' eccentricities, and the e\ resonance which per- 
turbs only the inner satellite's eccentricity. Subscripts 1 and 2 
represent Remus (inner) and Romulus (outer), respectively. 

The relevant resonance argu ments (j> for these three 3: 1 ec- 
centricity resonances are (e.g., |Murray & Dermott|1999 i 



4: 



e\e 2 



</» = 3A 2 -A 1 -2ro 2 , 
:3A 2 -A!-a3i-C32, 



> = 3A 2 -A,-2GJ I 



(10) 
(ID 



(12) 



where A is the mean longitude and G3 is the longitude of peri- 
center. Occupation of any of these resonances requires libra- 
tion of the considered resonance argument (exact resonance 
occurs when 4> = 0). Note that A = n, where n is the mean 
motion and is related to the semi-major axis. 



The resonance arguments given in Equations ( 10 1— ( 12 1 are 
listed in the order that these resonances are encountered due 
to tidal migration: first the e\ (at ai/a 2 ~ 0.481), then the eie 2 
(at ai/a 2 « 0.483), and lastly the e\ (at ai/a 2 ~ 0.485), where 
the resonance locations ai/a 2 must be adjusted depending 
on the exact starting values of a\ and a 2 . These resonances 
are not located at the same semi-major axes; differentiation 
of Equations (p~0[> — <fl2~|> shows that the various resonant argu- 
ments will librate atdifferent values of n\ and n 2 . Such "res- 
onance splitting" occurs because perturbations such as the ef- 
fect of primary 7 2 , and to a lesser degree (in this case, 2-3 
orders of magnitude smaller), secular perturbations, caus es Gi 
of the satellites to precess at different rates. In Section [6\4] 
we will discuss the capture of Remus and Romulus into any 
of these resonances. Next, we describe our methods regard- 
ing the implementation of tidal effects using direct N-body 
integrations. 

6.3. Methods 

Our methods and implementation for simulating a 3:1 res- 
onant passage due to tidal migration are as follows. We use 
an N-body integrator with a variable-timestep Bulirsch-Stoer 
algorithm from Mercury (Chambers 1999). We implement 
additional terms in the equations of motion due to the effects 
of tides on semi-major axis and eccentricity by fo llowing the 
numer ical methods described in Appendix A of Lee & Peale| 
(2002). Specifically, we used Equations (H} and (|6j) to model 
the tidal evolution in time of a and e. W e nave tested our im- 
plementation by reproducing results in |Lee & Peale ( 2002 



as well as matching the analytical expectations (Equations (|4 1 
and (|6]l) of semi-major axis drift and eccentricity evolution 
outside of resonance. 

Actual tidal timescales can be computationally prohibitive, 
and we incorporate a "speedup" factor to artificially increase 
the rate of tidal evolution in our simulations. Such speedup 
factors have also been numerically imple mented in previ- 
ous studies of tidal migration (i.e., Ferraz-M ello et al^2003| 
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Meyer & Wisdom|20 08 1 Zhang & Nim mo|2009) , where they 
adopted values up to 1000 and found that their results were 
not sensitive to the choice of speedup factor in the range 1- 
1000. In our implementation, we incorporate a speedup fac- 
tor by multiplying Equations Q and ([6]) by typical speedup 
factors of 100- 1000. In agreement witn previous studies, we 
find that our results are not sensitive to the choice of speedup 
factors up to 1000 for select test cases. 

We integrate the system for artificial durations of 1-10 Myr, 
which, because of the speedup factors, represent 1 Gyr of tidal 
evolution. Our figures show the tidal evolution timescale, 
not the artifical timescale used in the integrator. The 1 Gyr 
timescale is constrai ned by the lifetime of Sylvia's satellite 
system. Work by Vokrouhlicky et al. (2010) investigating the 
collisionally-born asteroid family related to Sylvia suggests 
that the family members (and hence the satellites) are at least 
10 8 years old. We can also estimate the lifetime of the satel- 
lites by considering how much time would pass before a col- 
lision between one of the satellites and another main belt as- 
teroid. We estimate this timescale to b e roughly 10 9 years 
( |Farinellaetal.|1 998 ; Bot tke et al.|2005| l by assuming that the 
smallest satellite has a diameter of 10.6 km, as suggested by 
our orbital fit analysis (Section [3J. Accordingly, we consider 
1 Gyr to be a reasonable time within which tidal evolution 
can have taken place, and we typically do not run simulations 
longer than 1 Gyr. 

For the default set of initial conditions in our simulations, 
the masses of all bodies and primary oblateness (J2) are taken 
from Table |4] Simulations are started with coplanar and 
nearly circular (e = 0.001) orbits. Angles for the argument of 
pericenter and mean anomaly are given a random value from 
0° to 360°. Initial semi-major axes for Remus and Romulus 
are 654 km and 1352.5 km, just inside the 3:1 resonance loca- 
tion. To simulate tidal evolution, we also need to adopt values 
for the Love number k and tidal dissipation factor Q. To cal- 
culate the Love numb er, we use the tidal mon olith model for 
these bodies (see e.g., Goldr eich & Sari||2009[ l. For all bod- 
ies, we assume Q = 100, a reasonable assumption for rocky 
monoliths. 

We also ran additional sets of simulations where we var- 
ied the initial eccentricities of both satellites (0.001-0.050), 
primary J2 (5-10% lower and higher than its best-fit value), 
satellite masses that spanned the range of adopted la con- 
fidence intervals (low-high, low-low, high-low, high-high 
combinations), and repeats of the nominal configuration for 
additional randomized initial angles of the argument of peri- 
center and mean anomaly. We did not specifically vary the 
tidal quantity k/Q, as the effect of varying k/Q is the same as 
varying the speedup factor since they both contribute linearly 
to d and e. We describe our results in the next subsection. 

6.4. Results: Evolutionary Pathways 

Here we describe the results stemming from our simula- 
tions of a 3:1 resonant passage between Remus and Romulus. 
These results suggest three evolutionary pathways: capture 
into resonance with no escape, temporary capture followed 
by escape, and no capture. We describe each of these evolu- 
tionary pathways in the following paragraphs. 

6.4. 1 . Capture with no escape 

This scenario, where resonant capture occurs with no es- 
cape during the 1 Gyr evolution, was typically observed for 
the e\ resonance. An example of such evolution is shown in 



Figure [7] Given that (a) the satellites are not currently ob- 
served in the 3:1 resonance and (b) these simulations show no 
escape from such resonant capture within a reasonable system 
lifetime of 1 Gyr, we conclude that this evolutionary pathway 
did not occur and we do not discuss it further. 
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Fig. 7. — Example of a simulation with resonant capture into 
e\ and no escape. Dots show results from numerical simula- 
tions for the libration angle <fi and eccentricity e as a func- 
tion of time. Subscripts 1 and 2 represent the inner and 
outer satellites, respectively. The simulation was started at 
ai/«2 = 654.0/1352.5 = 0.4835 and reached the e\ resonance 
at ai/a 2 = 655.2/1352.5 = 0.4844, roughly 8 Myr after the 
start of the simulation. The outer satellite's eccentricity will 
continue to increase due to the resonance effects, but the tidal 
damping effects will increase as the eccentricity grows, such 
that an equilibrium value for e may be reached. 



6.4.2. Temporary capture followed by escape 

In this event, resonant capture occurs and is followed by 
eventual escape due to growth of the resonant argument. This 
was a common outcome for each of the 3 types of resonances. 
Examples of such evolution are shown in Figure [8] Final ec- 
centricities at the end of our simulations ranged from their ini- 
tial values up to ^0.3. We are unable to place lower bounds on 
the final eccentricities because we cannot assess how long the 
satellites may have been captured in the resonance. If high ec- 
centricities resulted from temporary resonance capture, then 
sufficient eccentricity damping may have subsequently oc- 
curred to bring high post-resonance eccentricities to low ob- 
served eccentricities. 

We investigated whether such eccentricity damping was 
possible within a conservative timeframe of 1 Gyr. To do 
so, we integrated the a and e tidal expressions in Equations 
|4} and J6} and considered all Love number models (Equa- 
tions VWj, various post-resonance eccentricities up to 0.25, 
and tidal dissipation Q values (10-1000). We assumed that 
both satellites had densities equal to that of the primary (1.29 
g cm -3 ). 

For Remus, we find that eccentricity damping to observed 
values is only possible if we assume rubble pile Love num- 
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ber models (either k oc R or k oc l/R) for Remus (there is 
no restriction on the primary). If we make this assumption, 
damping to observed values is possible by adopting reason- 
able values of Q p =100 and Q s = 10-100. When we assume 
that Remus is monolithic, even when we adopt very favorable 
conditions for eccentricity damping tidal damping to its ob- 
served value is only possible if we assume a post-resonance 
eccentricity of ^0.032 or less. These calculations suggest that 
if its post-resonance eccentricity exceeded ^0.032, it is likely 
that Remus may have an interior structure more akin to a rub- 
ble pile aggregate than a monolithic body. 

For Romulus, damping to observed eccentricities is possi- 
ble only if the eccentricity was barely affected while in the 
resonance (as well as assuming favorable dissipations condi- 
tions: Q p = 1000 and Q s = 10). If the eccentricity reached 
even modest values (~ 0.023) we find that none of the Love 
number models and reasonable 2=10- 1000 values can damp 
eccentricities to even the highest possible observed eccentric- 
ity (0.01 1) allowed by our fit uncertainties. Therefore, if tem- 
porary capture in the 3:1 occurred, it must not have lasted 
long enough for the eccentricity of Romulus to reach values of 
~ 0.023. While such a scenario does not entirely rule out the 
e\ and e\e 2 resonances, it does seem to place bounds on the 
acceptable increase in eccentricity due to the 3:1 resonance. 

6.4.3. No capture 

In this case, all eccentricity-type resonances are encoun- 
tered and none result in capture. For orbits that are slowly 
converging toward each other, capture is possible depending 
on their initial eccentricities. When the pre-encounter eccen- 
tricity is below a critical eccentricity, capture is guaranteed. 
When the pre-encounter eccentricity exceeds a critical eccen- 
tricity, capture becomes a probabilistic event. For t he 3:1 res- 
onance, crit ical eccentricities can be estimated as (Murray & 
|Dermott|1999l > 



<?l,crit : 



32/i 



3 2 / 3 ^ + „ 

M 2 M 2 M 2 
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where subscripts p, 1, and 2 represent the primary, Remus, 
and Romulus, respectively. The fi and f 2 terms represent 
functions of Laplace coefficients b^ 2 (oi). They are (e.g., see 
|Murray & Dermott|1999| l 



1 



fi = -(-5j+4f-2aD+4jaD+a 2 D 2 )b (j Ua) 
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r 

-a, 
(16) 



where j = 3 for the 3:1 resonance, a = a\/a 2 is the ratio of 
semi-major axes, and D = d/da is the differential operator. 
The quantity -27 a/8 in Equation ( 16 1 is the indirect term for 
the case when M 2 > Mi . 

Using Equations ( p~3] > and (fl4) i, we calculate the critical ec- 
centricities to be ei.crit = 0.00864 and e2,crit = 0.00410. Even 

3 Qp — 1000 and Q s = 10. Inspection of Equation J&J shows that damping 
can be speeded up by making Q p / Q s as large as possiBle. 



for initial ej values that are low (e.g., osculating value of 
0.001), ei < ei.d-it will not always be satisfied because short- 
term eccentricity fluctuations due to primary J 2 will inflate ej 
excursions up to ^0.014 (assuming cii = 654 km in Equation 
([3]l of Section BJ. For Romulus, marooned farther from the 
primary such that J 2 effects are lessened, if the pre-encounter 
eccentricity is low enough then it is possible that the eccen- 
triciy will always remain less than the critical eccentricity. 
These analytical arguments are in agreement with the results 
from our numerical experiments. 

If we contemplate scenarios in which resonant capture 
never occurred in Sylvia's past, then we must adopt the critical 
eccentricities as lower limits on the past eccentricities of Re- 
mus and Romulus. Their past eccentricities cannot be lower 
than these limits because otherwise capture would have been 
guaranteed. We note that these lower limit eccentricities are 
lower than the nominal observed eccentricities (Table [4]), and 
hence this evolutionary pathway is a plausible scenario with- 
out requiring any significant modifications in eccentricity over 
time. 

To summarize these results, from these 3 pathways we find 
that both (a) temporary capture followed by escape and (b) 
no capture are plausible scenarios that occurred when Re- 
mus and Romulus encountered the 3:1 resonance. If pathway 
(a) occurred, our calculations of the necessary damping re- 
quired to bring post-resonance eccentricities to observed val- 
ues show this is possible for Remus but may be prohibitively 
long for Romulus, depending on its post-resonance eccentric- 
ity. Therefore it is unlikely that a substantial increase in the 
eccentricity of Romulus occurred, even if the system was tem- 
porarily captured in the e\ or e\e 2 resonances. If pathway (b) 
occurred, we can set lower limits on past eccentricities of both 
satellites to be equal to their critical eccentricities. 

7. CONCLUSIONS 

The goals of this study were to characterize Sylvia's current 
orbital configuration and masses as well as to illuminate the 
past orbital evolution of this system. Our work can be sum- 
marized as follows: 

(1) We reported new astrometric observations of Sylvia in 
201 1 that increased the number of existing epochs of astrom- 
etry by over 50%. These new observations extended the ex- 
isting baseline of observations to 7 years (for Remus) and to 
10 years (for Romulus). 

(2) We fit a fully dynamical 3-body model to the available 
astrometric data. This model simultaneously solved for or- 
bits of both satellites, individual masses, and the primary's 
oblateness (Table |4|. We found that the primary has a den- 
sity of 1.29±0.39 g cm" 3 and is oblate with a J 2 value in 
the range of 0.0985-0.1. Constraints on satellite radii can 
be obtained from the mass determinations by assuming that 
the satellites have a bulk density equal to that of the primary; 
we find ~4. 5-6.1 km for Remus and ~2.6-8.2 km for Ro- 
mulus. These ranges would have to be modified if the actual 
density of the primary or of the satellites was different from 
the nominal value assumed here. The orbits of the satellites 
are relatively circular. We find that the primary's spin pole is 
best fit when aligned to Romulus' orbital pole, and that the 
satellites' orbit poles are coplanar to within one degree. 

(3) We numerically investigated the short-term and long- 
term stability of the orbits of Sylvia's satellites. There are pe- 
riodic fluctuations in eccentricity for both satellites, most no- 
tably for the inner satellite Remus. We verified that these ec- 
centricity excursions are due to the effects of primary oblate- 
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Fig. 8. — Examples of simulations with temporary resonant capture followed by escape. Dots show results from numerical 
simulations for the libration angle <f> and eccentricity e as a function of time. Subscripts 1 and 2 represent the inner and outer 
satellites, respectively. Left-side plots: Results of a simulation with temporary capture into resonance, where both e\ and e 2 
increase. Initial conditions for eccentricity are e\ = 0.001 and e 2 = 0.04. Right-side plots: Results of a simulation with temporary 
capture into e\ resonance, where only e\ increases. Initial conditions for eccentricity are e\ = 0.007 and e 2 = 0.04. 



ness. From long-term integrations we found that the system 
is in a very stable configuration, in agreement with previous 
investigations. 

(4) We studied the past orbital evolution of Sylvia's satel- 
lites, including the most recent low-order MMR resonance 
crossing, which is the 3:1. We used direct N-body integra- 
tions with forced tidal migration to model such an encounter. 
To examine the case of resonant capture followed by escape, 
we calculate the tidal damping timescale to go from the post- 
encounter eccentricity to the observed value. Using available 
tidal models, we find that the damping timescale for Romu- 
lus can be prohibitively large if its post-resonance eccentricity 
exceeded ^0.023. This suggests that the system crossed the 
e\ and e\e 2 resonances without capture, or that it was not cap- 
tured in these resonances for a sufficient duration to substan- 
tially increase the eccentricity of Romulus. Similar timescale 
constraints from tidal damping also imply that Remus may 
have a rubble pile structure if its post-resonance eccentricity 
exceeded ^0.032. Alternatively, if no capture in any reso- 
nance occurred then we are able set lower limits on their past 
eccentricities {e x = 0.00864 and e 2 = 0.00410). 

The detailed characterization of Sylvia presented in this pa- 
per has allowed for analyses of its orbital evolution. Such 
studies of triple systems are important in order to understand 
their key physical properties, orbital architectures, and in- 



triguing evolutionary histories. 
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